Inference from Matrix Products: A Heuristic Spin Glass Algorithm 
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We present an algorithm for finding ground states of two dimensional spin glass systems based 
on ideas from matrix product states in quantum information theory. The algorithm works directly 
at zero temperature and defines an approximate "boundary Hamiltonian" whose accuracy depends 
on a parameter k. We test the algorithm against exact methods on random field and random bond 
Ising models, and we find that accurate results require a k which scales roughly polynomially with 
the system size. The algorithm also performs well when tested on small systems with arbitrary 
interactions, where no fast, exact algorithms exist. The time required is significantly less than 
Monte Carlo schemes. 



The problem of finding spin glass ground states is a 
major problem in statistical physics. Even for a two di- 
mensional Ising system, this problem is NP-hard if the 
bonds and fields are arbitrary[l]. For this reason, exact 
algorithms will often be unable to find the correct ground 
state in a reasonable time, and heuristic algorithms be- 
come necessary. 

Monte Carlo algorithms, including parallel 
tempering[2], extremal optimization [3], and others, 
present one important class of heuristic algorithms. 
Their performance can be tested against exact algo- 
rithms in two important cases, the case of planar Ising 
models where all the fields vanish but the bonds are 
arbitrary and the case of ferromagnetic Ising models with 
arbitrary magnetic fields, where fast exact algorithms 
are available. These problems are the random bond and 
random field problems, respectively. 

Unfortunately, these Monte Carlo algorithms can still 
take exponentially long time to find the true ground 
state even in these cases[4, 5]. In this paper, we 
present an approach to finding spin glass ground states 
inspired by applications of density matrix renormal- 
ization group (DMRG)[6] to two dimensional classical 
thermodynamics [7]. We call our approach the inference 
matrix product (IMP) algorithm. IMP, working directly 
at zero temperature, greatly reduces the time required 
compared to transfer matrix DMRG for these systems. 

We consider the following Hamiltonian: 



E = 2_^ JijCTiOj + hiOi, 



(1) 



<ij> 



where ai — ±1. The couplings Jij are between nearest 
neighbor spins on a two dimensional lattice (we consider 
square lattices in this paper, but the development is sim- 
ilar for other lattices) . We consider a system of length L 
and width W , with open boundary conditions, with the 
coordinate in the length direction labeled from 1 to L. 

We begin by presenting the algorithm for this problem. 
We then compare to exact algorithms on certain cases, 
and discuss the possible role of criticality in describing 
these systems. 



Inference Matrix Product (IMP) Algorithm — In this 
section, we present the IMP algorithm. We begin by 
some background, giving the solution of the problem in 
the case W ^ 1. In other words, we begin with a one 
dimensional chain. Let i?("^(i7i, ...jan) denote the energy 
of the first n spins for a subchain of length n < L. Let 
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CTl, 
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(2) 



denote the minimum energy of a subchain of length n < L 
with the last spin fixed to have value r. The subscript 
"^" denotes that the energy is over a subchain to the left 
of the spin we consider, and later we introduce energies 
E^; these energies are boundary Hamiltonians describ- 
ing the energy as a function of spins at the boundary of 
the chain. Then we have the recursion relations: 

EliTJ^h, (3) 

^(n+l)(^) = min,(^(r)(a) + Jn,n+1<JT + /i„+it) . 

These recursion relations can be solved in time 0{L) to 
find E^^t), which can be minimized over r to find the 
ground state energy of the chain. 

Similarly, we can find the ground state as follows: by 
minimizing El^'{aL), we can fix a^. Then, we minimize 
Ei^~^K<^L-i) + Jl-i,l<tl-iO'l to fix (Tl_i. Proceeding 
in this fashion moving back to the left along the chain, we 
find the ground state. The approach we have described 
is variously referred to as transfer matrix, dynamic pro- 
gramming, or belief propagation in different fields. 

To solve the problem with W > 1, we define 
Ei^'{ti, ...,Ty\/) to denote the minimum energy of the 
subsystem with length n and width W , with the last col- 
umn of spins fixed to be ri, ...,tw- We use coordinates 
(r, c) to denote the row and column of the spins, with 
r = 1, ...,L and c—1, ..., W. Then, 



sW(ri,...,TM/) 
i?(l'+')(ri,...,TH.) 



nima^^...,a„(E^J^\ai,...,aw) 



w 



(4) 






where the energy of a column of spins C'^"-' 



^w-i 



Err — i T . 



C(")(Ti,...,rv^) 

Unfortunately, the brute force solution of Eqs. (4) is 
not practical. There are 2^ possible spin configurations 
on a given column, and hence the time required scales as 
0(L22W'). 

The approach we follow is designed to solve this 
problem. Our approach, inspired by the idea of ma- 
trix product states, is to define an approximation to 

E^^Hti, ...,Tw), which we denote E]^ f.{Ti, ...,tw) and 
which we refer to as a matrix product energy. The quan- 
tity k is an integer, which we call the bond dimension. 
For larger values of k this approach becomes more accu- 
rate, but at the same time more computationally costly. 
This approximation will have the property that 



E]:^\{Ti,...,Tw)>E^:iHn,...,Tw). 



(5) 
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We define the energy E]^'^{ti^ ...^tw) by introducing 
auxiliary variables, called bond variables, on each verti- 
cal bond in the column of spins. We label these bond 
variables ai^i+i, for z = 0, ..., W, and each bond variables 
may assume k different values. Note that although we 
consider open boundary conditions for the interactions, 
we introduce bond variables ao,i and aw,w+i connecting 
above the top spin and below the bottom spin. This is 
done for notational convenience later. 

We now define an auxiliary problem which is a one di- 
mensional problem of the spins r and the bond variables 
a interacting along the column, with energy F^"' , and we 
define the energy E^ j.{ti, ..., tw) to be the minimum of 
pin) Q^Qj. ^YiQ bond variables a. Specifically, we set 



w 

pi^) ^ Y" F, 

i=l 



in) 



(Tj,Q;j_i^i,Q;j^i+l), 



(6) 



where the functions Fi are described below. Thus, each 
spin r is coupled to the bond variable above and below 
it. We set Jo,i — Jw,w+i — 0. 

Note that the energy e2:^{ti, ...,tw) from the first 
line of Eq. (4) can be described exactly by a matrix 
product energy with /c = 2 as follows. To see this, 
label the two states of the bond variable by ±1, just 
as we label the states of the spin by ±1. Then, set 
Fi (Tj,ai_i,i,Q;j,i+i) = /ii'ri + J(i_i_i)^(i^i)Q;j_i^iTj if Tj = 

ai^i+i and F/ ^(r^, ai_i,i, a^^i+i) = cx) otherwise. 

It is convenient for numerical purposes to "encode" 
the spins in the bond variables; that is, we will con- 
sider only states such that given a^^i+i, there is only 
one choice of r^ which gives a finite energy. Then, 
given a matrix product energy, E]^ i,{ti, ...,tw) with 
bond dimension k, we can find a matrix product 
energy E]^ ^i^{ti, ...,tw) which has bond dimension 
2k and which equals miners <jvi,(£'!l^fc(o'i, •••,o'vf) + 



Ej^i Jinj)Xn+i,j)CrjTj + C("+iHn, ■■■,Tw)) ■ We refer 
to this as "propagating" the energy to the right. 

As this propagation proceeds, k increases exponen- 
tially. Thus, an exact way to compute the ground state 
energy is to proceed as follows: initialize the energies 
Fj to give a matrix product energy with fc = 2, and 
propagate to the right, until the final matrix product 
energy has bond dimension k\y — 2^. Then, the task 
of finding the optimum spin configuration proceeds as 
follows. First, compute the optimum spin and bond 
variable configuration for the one dimensional problem 
defined by energy F^^>; since this is a one dimensional 
problem, the minimization can be readily done in a time 
0{k^W) by the techniques described above for one di- 
mensional chains. Then, fix the spins on column L to 
the values which minimize F^^\ and compute the opti- 
mum spin and bond variable configuration for the prob- 
lem F(^-i)-hX;j=i -^(L-i,»),(L,i)CT(L_i,i)tT(L,i). Proceeding 
in this fashion, we find the ground state. 

Of course, this approach is also exponentially costly. 
Therefore, we define below a removal of redundances 
and a "truncation" . We fix a maximum bond dimension 
kmax, and when the bond dimension of a matrix product 
energy F|!^\, exceeds kmax we truncate by finding a good 

" (n) 

approximation E^ j, . This truncation is achieved by 
choosing, for each bond, a subset of size kmax of the bond 
variables, and only allowing the bond variables to assume 
those particular choices. Then, we will have the property 
tnat lor k ^ i^mn^i 



E' 



in) 



> E 



(n) 



(7) 



The way we choose this subset is described below. Un- 
surprisingly, the best choice of the bond variables will 
depend on the optimum configuration of the spins to the 
right of those we have considered thus far, and at this 
point we have no knowledge of the best configuration of 
those spins. Therefore, after propagating E^^k to the 
right and truncating, we then propagate back E^^k to 
the left so that on the propagation back to the left we 



will use the E 



(n-l) 



which we have previously computed 



to decide the best truncation of E 

The IMP algorithm is defined as follows: 

1 . Initialize the energies F^ and set the bond dimen- 
sion for the first column, fci, equal to 2. Set n — I. 

2. Propagate to the right to compute -E^\ from 

" (n) 

El_i^, and set the bond dimension kn+i for the 
n + I-st column equal to 2A:„,. 

3. If kn+i > kmax, for somc kmax, rcmovc redundan- 
cies and truncate as described below. 

4. Increment n and if n < M^ goto step 2. 



5. Repeat steps 1 to 4, but this time propagate E_^\. 
from right to left. 

6. Compute the optimum spin configuration, working 
from left to right. 

Removing Redundancies — Before truncating, we re- 
move redundancies. In practice, this step does not have 
much effect on performance, and so on first reading the 
reader should skip down to the section on truncation. 
The bond variable q;o,1) leaving the top spin connects to 
only one spin, ri. Therefore, we can replace any matrix 
product energy in which this bond variable can assume 
k different variables by one in which the bond variable 
can assume only a single value, which we set to zero, by 
replacing Fi(ti, ao,i, ai,2) by F' with F'(ti,0, 0:1,2) = 
min^p j^i^i(ri, q;o,i, q;i_2)- We can do something similar 
for the bottom variable, aw,w+i, except since this vari- 
able encodes the state of tw , we must allow it to assume 
two different values. 

For all the bond variables in between, we can 
remove redundant states as follows. Consider 

a bond variable a^^i+i. If there are two states 
of this bond variable, h and l2, which have the 
property that i^j(Ti, a^-i^i, ^i) = Fi{Ti,ai-iA,l2) 
for all Qfi-i^i, then we replace -Fi+i by F' with 
F'{Ti^i,ai^i+i,ai+i^i+2) ~ ^i+i('''i+i, Qfi^i+i, ai+i,i+2) 
if Qfj^i+i ^ h and F'(ri+i, Zi, 0:^+1,^+2 = 

minf Fi+i (n+i, ^i, Q;j+i,i+2), -Fi+i (r^+i , I2, aj+i,i+2) j • 
Now ai^i+i is only allowed to assume fc — 1 dif- 
ferent values and no longer may assume the 
value l2- We can do a similar procedure if 

Fi+l{Ti+i, h, tti^i+i) = Fi+i{Ti+i,l2, Oti^i+l). 

Finally, if it turns out to be the case for a bond variable 
ai^i+i we find that Fj(tj, ofj-i^i, ^i) < i^i(ri, ai_i,i, ^2) 
for all Ui^i^i and all r^, and that i^j_|_i(rj_|_i, Zi,ai,i_i_i) < 
Fi+i{Ti+i,l2,at^i+i) for all tti^i+i and all n+i, then we 
can remove the choice l2- 

Truncation and Parameters — We now describe the 
truncation procedure. We describe the procedure in the 
case of the propagation back to the left (step 5), where we 
truncate a matrix product energy Ei^'^. To help guide 
the choice of truncation, we use the matrix product en- 
ergy Ei^^, . We will write F^ and F^ to refer to the 



Fur- 



auxiliary problems used to define E^ ^, and Ei^ j, . 
ther, we write Ji to refer to J{n-i,i),{n,i)i and we use al 
to refer to spins in column n and a\ to refer to spins in 
column n — 1, and use ai^i+i to refer to bond variables 
in F^ and /3i,i+i to refer to bond variables in F^. 

We begin by truncating bond variable aw,w+i- We 
consider the energy 



Et. 



F^ + F^ 



^J,a: 






(8) 



For each of the k different choices of aw,w+i, we min- 
imize this energy over the other bond variables a and 



/? and the spins. We then truncate by only keeping the 
kmax values of aw,w+i which lead to the lowest energy. 
We then truncate the bond variable aw~i,w by com- 
puting, for each of the k different choices of this bond 
variable, the minimum of the energy over all the other 
the bond variables and spins, and again keep only the 
kmax values which minimize the energy. We proceed in 
this way until we have truncated all the bond variables. 
Note that in, for example, the truncation of aw~i.Wi 
when we minimize the energy we only allow bond vari- 
ables aw,w+i to assume one of the kmax different states 
we kept on the first truncation step. 

To do this truncation requires computing the energy 
the minimum energy of a one dimensional system, Efot- 
This can be done by computing the minimum energy of 
a subchain consisting of the bottom h bonds, where the 
top two bonds are constrained to have values a and /3. 
Propagating this energy, just as in Eq. (3) gives the de- 
sired answer. The most naive algorithm would take a 
time 0{Wk'^) to do the truncation. However, we can 
take advantage of the simple form of the interaction be- 
tween the two chains, and do the entire truncation step 
in a time 0{Wk^), by propagating the energy upward for 
each of the four choices of the top two spins. 

In step (3), we perform the truncation similarly to here, 
except that we set Et^t = F^, without considering the 
interaction between two columns of spins. While this is 
much faster, taking only a time O(W^fc^), the propagation 
back in the other direction in step (5) is necessary in 
many cases to obtain accurate results. The total time for 
the algorithm then scales as 0{LWk^). 

Results on Random Bond and Random Field Ising 
Models — We compared the performance of this algo- 
rithm to exact solutions on the random bond and random 
field Ising models. In the random bond case, we chose 
each bond strength independently at random from a uni- 
form distribution between —1 and 1. All the magnetic 
fields were set equal to zero. While polynomial algo- 
rithms exist for planar Ising models with no fields, as in 
this case [8] , but we used the spin glass server to find exact 
solutions[9] for problems of size L = W = 5,10,20,40. 
In Fig. 1(a) we show the performance of the algorithm as 
a function of k for these sample sizes. As a rough esti- 
mate, the bond dimension for the random bond problem 
to obtain 50% accuracy scales as L^/^ ^ so the total time 
computational time required scales roughly as L^ . More 
numerical work is needed to truly ascertain the scaling 
of k with L, but certainly we do not need exponentially 
large values of L to obtain good accuracy. 

In the random bond case, we set all bonds ferromag- 
netic with unit strength, and chose the fields indepen- 
dently at random from a uniform distribution between 
—H and +H for some parameter H . For small values of 
H , we tend to produce only a single domain for small sys- 
tems, and so finding the ground state is simple (however, 
for any fixed H , the ground state will acquire domains for 



<; 0/ 





FIG. 1: Fraction of correctly identified ground states as a 
function of k for random bond (a) and random field (b) prob- 
lems (50 realizations in each case). The curve with fc = 5 is 
barely visible in a since perfect accuracy out of 50 samples 
was obtained with A; = 3. A small number (50) of samples 
was deliberately chosen to illustrate that for some individual 
samples increasing k can occasionally lead to a loss of accu- 
racy. 



may scale exponentially with L. However, for typical 
problems, generated randomly, a polynomial scaling of k 
is quite possible. This may be related to possible confor- 
mal invariance[16, 17] in the random bond problem, given 
the relationship between entanglement entropy and the 
scaling of k for one dimensional quantum systems[12-14]. 

It would be interesting to extend to three dimen- 
sional problems, by defining a boundary Hamiltonian on 
a plane, rather than a column. Truncating the bond 
variables would require minimizing the energy of a two 
dimensional problem, which can be done using the tech- 
niques described here. This idea of descending to one 
dimension lower is reminiscent of [18]. Another possi- 
ble extension would be to combine the algorithm with 
Monte Carlo methods, by using randomness in the choice 
of states to keep in the truncation, and doing several 
sweeps until the ground state is found. 
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by U. S. DOE Contract No. DE-AC52-06NA25396. 



large L[10]). For large values of i?, each spin aligns with 
the local field and again finding the ground state is sim- 
ple. The hardest cases seem to arise for moderate values 
of H . In Fig. 1(b), we show simulations of accuracy as a 
function of k, so systems sizes L — W = 10,20,40 with 
H = 5. This value of H was chosen to provide a system 
with large domains. Exact solutions were performed by 
an open source push-relabel algorithm[ll]. 

Finally we tested on systems with arbitrary fields and 
bonds. Exact solutions were available for small system 
sizes by brute force; for larger system sizes we can use 
the IMP algorithm for large k to test the algorithm for 
small k. The accuracy of the algorithm as a function of k 
is comparable in this case to the random bond and field 
problems, as far as we were able to test it. 

Discussion — We have presented an algorithm for find- 
ing ground states of two dimensional glassy systems. The 
sweep back left is important for accuracy. Working at 
zero temperature is important for speed; doing transfer 
matrix DMRG[7] or TEBD[15] for a system at non-zero 
temperature would require a time 0{k^) to do the initial 
sweep to the right, and would require a time O(fc^) in 
the sweep back left, compared to O(k^) and 0{k^) re- 
spectively for IMP. Note that the sweep back left to im- 
prove truncation was not considered in those algorithms 
since they were not applied to glassy systems. Running 
for k — 20 with L — W — 40 requires roughly 3 seconds 
of CPU time on a 2.0 Ghz G5, much faster than Monte 
Carlo. Given that for arbitrary fields and bonds the prob- 
lem is NP-hard, the required k for an arbitrary problem 
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